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Abstract. It has been known for a long time that the large experimental singlet-octet 
mass gap in the pseudoscalar meson mass spectrum originates from the anomaly of the 
axial vector current, i.e. from nonperturbative effects and the nontrivial topological 
structure of the QCD vacuum. In the Ncoiour oo limit of the theory, this con- 
nection elucidates in the famous Witten-Veneziano relation between the r;'-mass and 
the topological susceptibility of the quenched QCD vacuum. While lattice quantum 
chromodynamics (LQCD) has by now produced impressive high precision results on 
the flavour nonsinglet hadron spectrum, the determination of the pseudoscalar singlet 
mesons from direct correlator studies is markedly lagging behind, due to the compu- 
tational complexity in handling observables that include OZI-rule violating diagrams, 
like the rj' propagator. In this article we report on some recent progress in dealing with 
this numerical bottleneck problem. 



1 Introduction 

Long before the advent of QCD as the field theory of strong interactions, the 
TT-meson has been allocated the role of the Goldstone Boson in a near-chiral 
world of light quarks while its flavour singlet partner, the rj' meson, was thought 
to acquire its nearly baryonic mass of 960 MeV through ultraviolet quantum 
fluctuations that prevent the flavour singlet axial vector current 

jf\x)^q{xh,rq{x) (1) 

from being conserved, even in the massless theory. In this scenario the breaking 
of SU(3)l X SU(3)ii chiral symmetry occurs as a renormalization effect within a 
Wilson expansion of the operator that suffers operator mixing with the 

topological charge density 

Q{x) = igV32n^) trFixT'^F^^ix) (2) 

in the form 

d^jt;^°\x)^2NfQ{x) ■.= Q{x) . (3) 

The pseudoscalar operator on the right hand side of this ABJ anomaly [1] equa- 
tion is proportional to the number of active quark flavours, Nf and contains 
the gluonic field tensor, F, and its dual, F^ . Its Nf dependence is indicative 
of "disconnected" quark loop contributions as they appear in the socalled tri- 
angle diagrams of perturbation theory, see Fig. (1). Disconnected diagrams are 
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disconnected insertion of y_5 yji 




Fig. 1. Perturbative triangle diagram of flavour singlet axial current as disconnected 
diagram with 757^^ insertion. 

peculiar to matrix elements of flavour singlet operators and can be character- 
ized by intermediate states free of quarks, i.e. by quark-antiquark annihilation 
into gluons with subsequent pair creation processes. Their phenomenological 
role has been studied e.g. in various strong interaction fusion processes where 
light quark antiquark pairs annihilate strongly into charm (or bottom) quark 
antiquark pairs: 

u + u c + c or b + b . (4) 

These processes, as depicted in Fig. (2), are said to violate the OZI-rule [2]. 

Inspection of the renormalization procedure for this triangle diagram reveals 
a clash between quantum physics and symmetry requirements: the renormal- 
izazion program poses the alternative of either giving up gauge invariance or 
chiral symmetry. Accordingly, Eq. (3) reflects nonconservation of chiral symme- 
try, for the sake of preserving gauge invariance, in form of the AB J-anomaly (for 
an excellent review on the subject, see the Schladming lectures of Crewther [3]). 

The lesson to be learnt is that pseudoscalar flavour singlet mesons are deeply 
affected by and thus indicative for the topological structure of the QCD vac- 
uum. Needless to say that the time averaged vacuum expectation value of Q = 
Q/{2Nf) itself is zero, (0|(3(a;)|0) = 0, by translation invariance and parity con- 
servation of strong interactions^. As a consequence, the net topological charge 
of the QCD vacuum state is zero. The key quantity describing topological fluc- 
tuations of the QCD vacuum is the topological susceptibility, X: which is defined 

by 

X = J d^x(0\ g(x)Q(0) |0) . (5) 

In the so-called 't Hooft limit of the theory, when the number of colours, Nc 
is taken to infinity {Nc — > cxd at fixed value for the number of active flavours, 
Nf, in the weak gauge coupling limit ~ 1/Nc), Witten and Veneziano were 

^ There is no spontaneous breaking of parity in strong interaction physics. 
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Fig. 2. OZI rule violating quark fusion process with disconnected quark loops and 75 
insertions. 

able to derive a simple relation for the r/ mass in the chiral limit [4] 



where F^r denotes the pion decay constant. Within the Witten-Veneziano model 
assumption, that the real world is well described by the above limit relation, 
LQCD can predict the 77' mass by computing the topological susceptibility in 
quenched QCD. 

The lattice implications of Eq. (6) have been discussed in quite some detail 
in a recent paper of GIUSTI et al [5] . In this contribution we will focus on those 
problems that one encounters in the direct approach to extract physics from 
fermionic two-point correlators in the singlet sector. 

2 Prolegomena: pseudoscalars in LQCD 

In this section we shall consider a symmetric world with three flavours to enter 
the discussion on the computational tasks to be met. 

2.1 Flavour octet sector 

In the flavour octet sector of QCD (with mass degenerate u,d,s quarks), the 
nonconservation of the axial vector current is determined by the PCAC relation: 



which is much simpler in structure than its flavour singlet analogue, Eq. 3. In 
LQCD, the computation of the 7r-meson proceeds by hitting the vacuum state 
|0) with the quark bilinear pseudoscalar operator as given on the right hand side 
of Eq. (7). This creates both pions and excited 7r-like states. Let us consider such 
an operator carrying negative charge 



mi, — lim 




^Iquenched ; 



(6) 



(7) 



u{x)^zd{x)\Q) = \tt )x + excited states 



(8) 
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Fig. 3. The connected propagator Gp^{x,x') for pseudoscalar meson states with neg- 
ative charge, with a quark d and an antiquark u running from source x to sink x' in 
the unquenched QCD background field, with 75 insertions. 

with up and down quark operators, u and d respectively. The propagator in 
the pseudoscalar meson channel from x to x' , comprising the 7r~ -meson and its 
excited friends, would hence read hke this: 

G^p^^{x,x') = {0\d{x'h5u{x')u{xh5d{x)\0) . (9) 

This vacuum expectation value is most easily computed by Wick contracting 
the two valence quark operators with matching flavour, and the result can be 
simply expressed in terms of the two quark propagators, P„ = M~^{x,x') and 
Pd ~ Mj^{x,x') as follows: 

Gg^(x,x') = TT[P:{x,x')Pd{x,x')] , (10) 

where the trace is to be taken both in spin and colour space. In LQCD, the quark 
propagators (from source xq to sink x') can be computed in the background field 
of the QCD ground state (the vacuum) by solving the lattice Dirac equations 

MDix')Pix,x')^6{x-xo) , (11) 

where Md is the discretized form of the Dirac operator including the gluon 
interaction and the r.h.s. denotes a point source located at lattice site xq, see 
Fig. (3). This is done by applying modern iterative solvers, like BiCGstab [6]. 
Note that 

• Al£i{U} depends on the giuonic field configuration (i.e. the QCD vacuum 
configuration {U} that has been produced independently by a suitable sam- 
pling process such as Hybrid Monte Carlo [7,8]) 

• the efficiency of the solver is governed by the (fluctuating) condition number 
of the large sparse matrix M]j{U} [6], [9]. 

The simulation of the QCD vacuum on the lattice being a nonperturbative com- 
putation, the giuonic lattice vacuum configurations comprise the entirety of vac- 
uum polarization effects, i.e. any conceivable fluctation from creation and anni- 
hilation of quarks and antiquarks (sea quark effects). So in Fig. (3) all vacuum 
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Fig. 4. Gluon propagators including quark loop self energy diagrams. 

polarization diagrams as depicted in Fig. (4) are understood to be included. We 

note in passing, that the numerical computation of the complete inverse of the 
lattice Dirac operator is prohibitively expensive: its computational effort grows 
with lattice volume, which in todays simulations is anywhere between 16* and 
64* space-time lattice sites! So, in practical calculations of the 7r-meson, one re- 
stricts the source, xq, to just a few space points in a particular time slice, say in 
t = 0. The task for calculation is then the pseudoscalar two-point function: 

Gg^(a^o,a:') = (Tr[P„*(xo,x')Pd(xo,a;')]) , (12) 

where the brackets denote the average over an ensemble of some hundred gauge 
field vacuum states. Note that it is enough in practice to compute but a few 
columns of the entire inverse Dirac matrix per configuration for achieving suffi- 
cient accuracy of the estimate, Eq. (12). Once this pseudoscalar two-point func- 
tion has been calculated the pion contribution is obtained by removing the ex- 
cited states from the pseudoscalar channel. For this purpose, one first extracts 
the zero momentum contributions by summing over all (spatial) x' in time slice 
t': 

4's(* = 0'*') •= {Y.^4P:{xo,x')Pd{xo,x')]) . (13) 

x' 

This latter two-point function at large enough time separations t' between 
source and sink will be dominated by the ground-state of the channel, i.e. the 
pion. Prior to that, it contains a superposition of ground and excited states: 

C^^iO, t') = ^ c, exp{-Ei t') Co exp(-m^i') . (14) 

i 

It is obvious that the strict localization of the source at the very lattice site 
030 overly induces excited state contaminations; hence a suitable smearing of 
the source around this site will deemphasize such pollutions and will help to 
achieve precocious ground state dominance. This is welcome because the two- 
point correlator decreases exponentially and is prone to suffer large statistical 
fluctuations at large values of t' . 

The numerical task in exploiting Eq. (14) is then to optimize the smearing 
such as to achieve an early onset of a plateau: the latter is defined as the t'-range 



6 K. Schilling et al. 

over which Gpg(0,t') decays with a single exponential. This is tantamount to 
requiring the 'local pseudoscalar mass', mpg(t'), to remain constant: 

mpg{t') := -dtdn[c'^^^{0, t')] = const. = . (15) 

Smearing of sources or sinks is achieved by a diffusive iteration [8] process: 
think of replacing pointlike sources or sinks by 

5(xo) = (/.(°)(xo)^ 0(^+1) (x), (16) 

where ^^^^ {x) is peaked around xq and is constructed by the recursion 



0(^+1) (x) 



l + 6a 



(17) 



The index 'p.t.' stands for 'parallel transported', and the sum extends over the six 

spatial neighbours of x. As a result the entire scries of smeared 'wave functions' 
0(*\i = 1,2,---,A^ — lis gauge covariant under local gauge transformations. 
Quark operators are computed after N = 25 such smearing steps, with the value 
a = 4.0. The smearing procedure is applied to meson sources as well as to sinks 
for connected and disconnected diagrams. 

2.2 Flavour singlet sector 

In the case of SU{3) flavour symmetry the appropriately normalized flavour 
singlet pseudoscalar operator reads 

0^"\x) = ^Y.'l^{x)^^1^{x) . (18) 

Here the sum extends now equally over all (mass degenerate u, d, and s) 
quark flavours. The flavour singlet pseudoscalar meson propagator is estimated 
from the average 

G^^^^{x,x') = {0^°^Hx')0^°Hx)). (19) 

In this case the Wick contraction now induces both connected and disconnected 
(two-loop) diagrams, D, 

D{x,x') = ^tr[O(0)t(^')]tr[o(0)(a;)]^ , (20) 

such that finally 

Gf^{x, x') = Gf^{x, x') - NfD{x, x') , (21) 

as depicted in Fig. (5). Note that at this stage the trace operations in Eq. (20) 
refer to spin and colour space only and the Wick contraction procedure leads to 
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the weight Nf = 3 of the disconnected relative to the connected contribution, 

the latter being identical to the expression given in Eq. (13). The relative minus 
sign is due to Fermi statistics of the quark operators. 

We note in passing, that connected and disconnected correlators come along 
with intrinsic positive sign such that the above minus sign leads to a steepening 
of C(o)(i') w.r.t. C(s)(i'), iff D shows a soft drop in f. 




Fig. 5. The connected and OZI-rxiIc violating ('disconnected') contributions to the rj' 
propagator Gpg from Eq. (19). Note that tliese diagrams keep only track of the valence 
quark lines: all gluon exchanges due to the vacuum background field have been omitted 
for clarity of the graph. It goes without saying that in unquenched QCD, the gluon 
propagators implicitly induce sea quark loops as in Fig. (4). 

Let us mention that the terminology "two-loop correlator" might be mis- 
leading: it refers to the fact that we are dealing here with two (valence quark) 
loop insertions of the 75 operator. The unquenched QCD vacuum includes of 
course a host of sea quark loops that are im,plicitely taken into account when 
sampling unquenchced QCD vacua in a LQCD simulation. They would enter in 
form of vacuum polarization effects on the very gluon propagators as depicted 
in Fig. (4). 

It is illuminating to consider the mechanism for the origin of the singlet 
octet mass gap in a model to real QCD: let us start in the Nc = 00 world, 
where the Ua{^) symmetry holds and flavour octet and singlet mesons are de- 
generate, mo = mg. This degeneracy is only broken by the higher quark loop 
contributions which are suppressed by inverse powers in Nc- It is hence straight- 
forward to approximate the r]' propagator, starting out from this scenario, as 
an explicit multiloop expansion of quark antiquark annihilation diagrams in the 
pure gluonic background field of quenched QCD, as illustrated in Fig. (6). The 




Fig. 6. Pictorial receipe for setting up the 77's propagator in quenched QCD as an 
explicit multiloop expansion. 
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actual calculation can be greatly simplified by introducing an effective loop-loop 
coupling, /iQ. In momentum space the ?7'-propagator then takes the form of a 
geometric series: 



1 



+ ml 



-M5 



+ ml 



(22) 



The series is readily summed up with the result: 

p 1 

p' + {ml+i,l)' 



(23) 



which reveals that the effective loop-loop coupling is nothing but the mass gap 
between singlet and octet pscudoscalar masses: 



A 2 2 2 2 

Am =: ruQ — rUg = /Xq, 



(24) 



Wc emphasize that in this cfFcctivc model it is only by summation of the 
entire series to all orders in I that we arrive at a simple pole expression for 
the singlet propagator in momentum space and hence at a simple exponential 
falloff in time of the two-point correlator. Note that this siuriniation happens 
automatically in an unquenched QCD situation where all multiquarkloop effects 
are already implicitly included in the very gauge field vacuum configuration! 



3 The real world with nn ss mixing 

So far we have been leading a rather academic discussion in an SU(3) symmetric 
world. We would like next to describe the problem in a more realistic setting. 
This opens up a pandora box of four parameters for decay matrix elements plus 
the mixing angle among isosinglet particle states. Two degrees of freedom in the 
decay constants can be characterized as mixing angles [10]. 

3.1 Phenomenological approach: alignment hypothesis 

In order to ease the interpretation of flavour breaking effects in the isosinglet 
sector, Feldmann et al [11,12] have proposed some time ago an intuitive scheme 
to describe the rj — rj' phenomenology by use of a quark state Fock state repre- 
sentation in conjunction with an alignment assumption. The latter reduces the 
number of parameters to two couplings plus a single mixing angle of the particle 
states in that Fock space. Let us briefly review their scenario: they proceed from 
the two isosinglet states 



\r]„) =^n\uu + dd)/V2 + 
\r]s) = ^s\ss) + ... , 



(25) 
(26) 
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where Vn.s denote light cone wave functions ^. The key assumption of their 
phcnomenological approach is their hypothesis of alignment of the current matrix 
elements: 

(0|9i757;*%l%) = ^ijfPn = n,s . (27) 

Due to the QCD interactions the physical r] and r]' states will turn out to be 
mixtures of these Fock states: 

lt)= U{4>) { ^J V= f " { ) ■ (28) 

\i) J \\Vs) J \svn-4> cos(j) J \\t]s) J 

As a consequence of the alignment hypothesis, Eq. (27), the four decay con- 
stants, fp 

{0\jU\P} = (OlQa^^^P) = fpP^^ i = n,s; P = v,v' (29) 
can be expressed in terms of and the very rotation angle ^ [11]: 



By feeding the anomalous PC AC relation 

d^Qil^l^li = ^ruiqasqi + Q (31) 
into the four decay matrix elements 

{0\d%,\v) = r.ml {0\d%,\v') = f;,ml,; {i = n,s), (32) 
the authors of ref. [11-13] arrive at the mass matrix mixing relation 

with a single mixing angle that is in accord with Fock state mixing. In Eq. (33) 
the abbreviations stand for the explicit chiral symmetry breaking terms from 
nonvanishing quark masses, mi: 

lAi ■= -j^{0\Qa5qi\Vi); i = n,s . (34) 

In 

For consistency of the approach, they have to postulate the symmetry of the mass 
matrix, which sets a constraint among the matrix elements (0|Q|77i); i = n,s . 

The l.h.s. of Eq. (33) being a continuum relation, great care must be exercised 
in renormalizing the operators. Since the underlying anomalous PCAC relation, 

^ There dots in Eq. (25) represent higher Fock states, including \gg) and components 
from the heavy quark sectors, such as cc) . 
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Eq. (31) is actually to be seen as a Ward-Takahashi identity from chiral rota- 
tions, the standard Wilson lattice fcrmions provide a problematical scheme for 
regularizing the matrix entries in Eq. (33), with the Wilson lattice fermions lack- 
ing chiral symmetry at finite lattice spacing [14]. This complex situation would 
be largely improved when dealing with Ginsparg- Wilson fermions in unquenched 
QCD, whose lattice chiral symmetry prevents undesired divergencies from op- 
erator mixings [5]. On the lattice, the mass^ entries to the matrix Eq. (33) are 
best accessible from propagator studies in the singlet channels [15-18]. 

We would like to explain next, that the idea of Fock state mixing of eigen- 
states can in principle be pursued very naturally from a variational approach on 
the lattice, once we have fully unquenched QCD vacuum configurations in the 
future, without making reference to an alignment hypothesis. 



3.2 Mixing singlet mesons on the lattice 

On the lattice, we can prepare the the Fock states |r/„) and \r]s) by hitting the 
vacuum at some (smeared) source location, say x = (0,0,0,0), with the bilinear 
isosinglet operators, 

Oi(0)|0) =gi (0)759^(0)1 0), i = n,s (35) 

and by waiting long enough to see the ground states emerge. 

In the sense of a variational method, one can next work with a superposition 
of these operators to achieve earlier exponential decay of the correlator and/or 
lowering of the ground state, the ry-meson: 

Oa = cos aOn — sin a Og ■ (36) 

Suppose this variational ansatz delivers an immediate plateau for some value 
of the variational parameter, a = a*: in that case we have actually hit an 
eigenmode of the interacting system and the Fock state mixing angle is (j) = a*\ 
Since one expects little mixing with higher Fock states, one could then recover 
the 7]' ground state from the orthogonal combination 

= sin aOn + cos a Os ■ (37) 

In general the state mixing information is encoded more deeply in the corre- 
lation matrix (still in operator quark content basis) 

rM_fCnn{t)-2Dnn{t) V^Dns{t) \ , . 

Note that the weight factors 2 and \/2 in this relation originate from the combi- 
natorics when Wick contracting the degenerate u (with u) and d (with d) quark 
operators contained On 

qa^qn ■= -^(w75M + djsd) . (39) 
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If plateau formation does not occur immediatley, at tp = 0, one needs to 

follow the time evolution up to a later time step when the higher excitations 
have died out, and ascertain the plateau condition within the time interval tp —^ 
tp + At in form of a generalized eigenvalue problem [19] : 

^r. A ,\ ( + cosa\ , . , / + cosQ!\ 

C{tp^At)[^_.^^^^ =eM-m,At) Citp) ^.^^j . (40) 

The problem can be recast into a symmetric form 

C-^/\tp)C{tp + At)C-^/'^{tp) ■ ria. = exp(-m^Z\i) • J7„ (41) 
with the definition 

.„(y = c«/=(y(r32:) m 

According to the symmetric plateau condition, Eq. (41), the direction of the 

ground state vector, rjaitp), remains fixed under further time displacements At 
through the transfer operator, r(Z\i) = C-^/^{tp)Citp + At)C-^^^itp). Hence the 
rj mixing angle in the quark state basis is given by the direction of eigenvector 
T]a(tp) at the onset of the plateau and might differ appreciably from the varia- 
tional angle a. The problem Eq. (41) being symmetric, the 77' is retrieved as the 
second, perpendicular eigenstate of T(Z\t). 

It might be illuminating to put Eq. (41) into perspective with the previous 
plateau condition, Eq. (15): by Taylor expanding of C{tp + At) to lowest order 
in At one can readily recast Eq. (41) into a local condition, with a 'logarithmic 
derivative' matrix construct, C~^/'^{tp)C' {tp)C~^/^{tp): 

C-^/\tp)C'{tp)C-^/\tp) ■ noc{tp) = -mr, ■ Vccitp) . (43) 

In principle, the phenomenological alignment postulate of Eq. (27) can be 
tested on the lattice by computing another, nonsymmetric set of correlators 

4(t) = (0|S'^gi(t)757/.%(<)9j(0)75'?,(0)|0) , (44) 

which allows to extract the various decay matrix elements. Let us emphasize 
in concluding this section that the scenario presented here relates to fully un- 
quenched three flavour QCD; according to the discussion in 2.2 the mass plateau 
conditions can only be anticipated if the quark contents of the bilinear operators 
in Eq. (35) are equally present in the quark sea ^. 

4 The three computational bottlenecks 

In this section we shall focus in more detail on the numerical obstacles that one 

encounters when carrying out the above program. These difficulties arc compu- 
tationally much more severe than in nonsinglet spectroscopy and hence require 
the development of advanced algorithms and numerical techniques: 

^ see the discussion on the partially quenched approach in section 5.2 
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1. Present limitations in sea quark flavours. It was only in the sec- 
ond half of the nineties that LQCD had developed the means to tackle first 
semi-realistic large scale simulations of QCD beyond the valence (or 'quenched') 
approximation. The reason is that Wilson-like discretizations of the Dirac op- 
erator are mandatory once we wish to deal with strong interaction situations 
sensitive to flavour symmetry. But such Wilson-like lattice fermions need ex- 
tremely large supercomputer resources when it comes to actually simulate the 
effects of dynamical sea quarks, i.e. to sample unquenched CCD vacuum con- 
figurations. The hybrid Monte Carlo algorithm (HMC) is a nonlocal sampling 
technique which was developed and improved ever since the late eighties [7]; it 
continues to be the workhorse in all of the major QCD simulation projects with 
Wilson-like fermions [8,20,16,21]. For technical reasons, however, HMC has the 
shortcoming of being limited to even numbers of sea quark flavours. As a con- 
sequence, all the above projects neglected dynamical s-quarks in working with 
mass degenerate u and d sea quarks only, i.e. considered an SU(2) symmetric 
sea of quarks only. It is to be hoped that the next generation of large-scale QCD 
simulation with Teracomputers will overcome this limitation with new sampling 
algorithms. 

2. Trace computations. We have seen that the peculiarity of flavour sin- 
glet objects lies in the annihilation of the valence quark lines into the flavour 
blind gluonic soup which implies the need to compute disconnected diagrams. 
But inspection of Eq. (20) shows that the cc'-summation from momentum zero 
projection requires the evaluation of the trace of the inverse Dirac operator. This 
task is definitely beyond the reach of modern linear equation solvers for matrices 
of rank 10^ and more. Traditionally stochastic estimator techniques (SET) have 
been the popular workaround [22,23]. 

3. Noisy signals. The singlet spectrum being inherently determined by the 

physics of the quantiim flTictuations of the QCD vacuum, the two-loop correla- 
tors are bound to suffer from a serious noise level. This calls for the use of noise 
reduction methods, in order to circumvent the need for overly large ensembles 
of vacuum fi,eld configurations. With this motivation low mode expansion meth- 
ods have been proposed recently with considerable success to replace stochastic 
estimator techniques [15,24]. 

4.1 Stochastic estimate of fermion loops 

The basic idea of stochastic estimation of the momcntunm zero projected loop 
operator in Eq. (20), X^^, tr(75M~^(a;, t)) , is very simple: instead of attempting 
to solve Eq. (11) V on delta-like sources, one introduces socalled stochastic 
volume sources which are completely delocalized, on the entire lattice volume V: 

C=(a,6,---,ey)- (45) 

The components Xi are chosen to be real random numbers sampled from a normal 
distribution N{xi) = l/\/27rexp(— a;|/2) or socalled Z2 noise which is nothing 
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but equally distributed random numbers ±1. One solves the Dirac lattice equa- 
tion on an ensemble of Ng such random source vectors, {^"}, for each individual 
gauge configuration: 

C = J^M^'C ■ (46) 
and computes the ensemble average of the scalar products 

V 

(^,c))^ := EE^rcr ^tr[75m-i] . (4?) 

^ a— 1 i—1 

One might phrase this procedure as an all-at-one-stroke method, at the expense 

of nonclosed loops sneaking in, in form of contributions from j^M^^ with 

L J ij 

i ^ j. These undesired contributions are suppressed, however, due to the 'com- 
ponentwise orthonormality' relation 

a=l 

Note that Eq. (47) readily allows for restricting the grand trace operation over 
the entire space-time lattice, TR, to any particular time slice, simply by tnm- 
cating the i-summation to the appropriate three-dimensional subvolume, such 
as to yield the estimator for the momentum zero expression on time slice t, 

In principle, these considerations hold for all kinds of insertions, not just 75 
like in our example. In any case the 'sneakers' are suppressed as 0{N~^) terms, 
but with a strength that depends on the insertion. The practical advantage of 
this stochastic procedure is, that in the instance of 75 insertions one needs to 
compute only some few hundred solutions to Eq. (46) which is much less than 
0{V)\ But we should keep in mind that this is achieved only at the expense of 
injecting additional noise into the problem^. 



4.2 Low eigenmode approximation to fermion loops 

We start with the eigenvalue problem for the socalled Hermitian Dirac operator, 

Q5 ■■= I5MD , (49) 

* We would like to point out that some authors [25] have refrained from using noisy 
sources altogether by using a single 'solid' wall source with = 1 and relying on 

gauge symmetry to suppress the unwanted nonclosed loops, j^Mj^^ , withi / j. 

L J ij 

Actually this argument relies on Elitzur's theorem according to which the vacuum 
expectation vahic (O) vanishes for any such non gauge invariant operator O! How- 
ever, if the sample of vacuum gauge fields is limited to some few hundred independent 
configurations only this procedure offers too little control over these systematic er- 
rors; nevertheless one could refine this variant by applying the 'solid' wall source idea 
on an ensemble of stochastic gauge copies of each individual vacuum configuration. 
To our knowledge this has never been attempted. 
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which reads 



(50) 



The quark loop with 75 insertion at time sHce t is simply expressed in terms of 
this eigensystem: 

The disconnected two-loop correlator from Fig. (5) then has the form 



io i 



where f represents the time separation between source and sink and we have 
exploited translational invariancc by summing over all time locations of the 
source, to- This summation is for free, once the eigenfunctions are available^. 
The determination of the entire spectrum being of course prohibitively expensive, 
wc proceed by arranging the spectrum in ascending order of |Ai| attempting to 
saturate the eigenmode expansion, Eq. (52), with its lowest modes. The intuitive 
reasoning behind such a strategy is, that the underlying physics is expected to be 
encoded in the infrared rather than ultraviolet modes. After all, it is well known 
that for Wilson fermions 15/16 of the spectrum is related to the unphysical 
doublers that become frozen to order 0{a~^) as the lattice spacing a is going to 
zero. 




100 150 200 
Index of Eigenvalue 



300 



Fig. 7. Cumulative spectral distribution of the lowest eigenmodes at the smallest (lower 
curve) and largest (upper curve) sea quark masses of the SESAM project [24]. 

^ Wc mention that there might be other opportunities for spectral methods in LQCD 
with their potential to deal with all-to-all propagators: like the long standing string 
breaking problem or heavy-light meson scattering[26]. 
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Fig. (7) offers an impression on the cumulative distribution of the lowest 

modes of the Dirac operator in the situation of the SESAM simulation, at their 
smallest and largest sea quark masses. Note that the masses of valence quarks (in- 
serted into the Dirac operator) and sea quarks (inserted into the hybrid Monte 
Carlo sampling of gauge field configurations) arc chosen to coincide; Fig. (7) 
therefore reflects the feedback mechanism from dynamical effects onto the spec- 
tral representation as quark masses are diminished. Fig. (7) suggests that the 
important spectral activity is mostly within the lowest 50 modes. 
But how many modes do suffice for saturating Eq. (52)? 

0.016 I • • • • • 1 




50 100 150 200 250 300 
# of Eigenpairs 



Fig. 8. Pion correlator in spectral approximation, plotted versus spectral cutoff, Ac. 
Shown is a sequence of trajectories, as obtained on fixed values of t: the curves represent 
C^*) (t = 1) , C*^*' {t = 2), C**^ {t = 3) etc , from top to bottom [24] . In this representation 
the justification of truncation is signaled by a flat Ac-dependency. 

For the saturation of the two point correlators, one would anticipate the 
need for larger values for the cutoff, A'^, as one diminishes the time separations 
towards a few lattice spacings. This expectation is indeed confirmed by the 
numerical results for the pion corrrelator as displayed in Fig. (8): The diagram 
shows a sequence of curves which - in descendig order - refer to C{t = l),C{t = 
2), • • • , C{t = 16), which are plotted versus the spectral cutoff. For small values 
of t the cmves keep rising with Ac, while asymptotic saturation is found to be 
better at larger t-values. 

In Fig. (9), we plotted the pion correlator from a truncated eigenmode ap- 
proximation (TEA) with cutoff Ac = 300 and from standard iterative solvers. 
There remain noticeable differences in the two curves for this value of the cutoff. 
Therefore, one would refrain from using TEA for connected hadron propaga- 
tors unless really needed. The situation w.r.t. A^ turns out to be much more 
favourable with the disconnected piece, D{t), as shown in Fig. (10). 
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Fig. 9. Quality of the truncated eigenmode approximation at the smalles SESAM 

sea quark mass for the pion correlator, with a spectral cutoff, Ac = 300. Satisfactory 
agreement with the standard result from using iterative solvers is only reached at 
t ~ 10. 



How well the idea of saturation actually works out quantitatively is illus- 
trated in Fig. (11), where we compare the two-loop correlator in the truncated 
eigenmode approximation (TEA) against the results from a previous comprehen- 
sive study with the estimation by an ensemble of stochastic sources [23]. We find 
that the contributions from the lowest 300 eigenmodes agree remarkably well 
with the SET computations which are based here on 400 source vectors, on a 
16 X 32 lattice in the sea quark mass range covered used by the SESAM project! 
The net computational effort behind the two sets of data points in Fig. (11) is 
about equal, yet the TEA data among themselvesa appear to be less fluctuating. 

This positive message from Fig. (11) bears substantial promise for future 
lattice studies with lighter quark masses, as SET is doomed to degrade when 
penetrating more deeply into the more chiral regime while the truncated eigen- 
m,ode approxirn,ation (TEA) will thrive: (a) the lowest modes will dominate all 
the more as |Am,in| approaches zero; (b) modern Arnoldi eigensolvers have been 
made very efficient in present day LQCD simulations [27] and will not deterio- 
rate for lighter quarks as long as one sticks to fixed lattice volumes. Conversely, 
iterative solvers on stochastic sources will suffer substantial losses in convergence 
rate, once the condition number of the Dirac matrix, c := \Xmax/^min\ starts 
exploding. 

4.3 Navigating between Scylla and Charybdis 

There is another lesson from Fig. (11): we find the noise level in D(t) to be rather 
independent of t, i.e. of signal size. This is of course in line with the fact that we 
are hunting for vacuum fluctuations proper. It therefore appears impracticable 
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to achieve ground state projection in the very standard fashion by increasing t: 

the 77' correlator being the imbalance between connected and two-loop signals 
one is to keep track of a signal in form of a rapidly diminishing difference between 
two quantities, with one being at fixed noise level. Given this situation we are 
challenged by complying to two conflicting needs: on one hand large t values 
are needed to deflate excited state contaminations and on the other hand the 
struggle for an acceptable signal-to-noise ratio impedes working in the large-f 
regime. Our strategy to deal with this conundrum of the 77' correlator, Eq. (21) 
is to contrive a procedure for circumventing this conflicting requirements. The 
receipe is the following: 

• suppress the excited states from Cpg^At) by using an exponential fit from 
the window of its mass plateau; 

• deemphasize the short range noise level inside D by means of the truncated 
eigenmode expansion, Eq. (52). 

This method could also be characterized through an intermediate 'synthetic data 
stage' since the actual ground state analysis for the singlet correlator is organized 
as a sequence of three steps: 

(a) extract an exponential fit to the lattice data for E{At) = Cpg{At), 

(b) combine this fit curve, E{At), with the lattice data for D{At) into the 
'synthetic' data set for C^pl{At) := E{At) - NfD{At), in the spirit of Eq. (21), 

5.5e-14 
5e-14 
4.5e-14 
.0 4e-14 
I 3.5e-14 
^ 3e-14 
-? 2.5e-14 
5 2e-14 
1.5e-14 
1e-14 
5e-15 

50 100 150 200 250 300 
# of Eigenvalues 

Fig. 10. Quality of low eigenmode approximation as function of time separation 
for the OZI rule violating term, D(t). The curves represent (from top to bottom) 
C7(1),C(2),.--,C7(16) [24]. 
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Fig. 11. The two-loop correlator, D{At), as obtained from SET [23] (circles) and TEA 
(triangles) in comparison. Local sources and sinks have been used and Ns = 300 and 
Ac = 300 were chosen [24]. 

(c) establish a mass plateau in Cpg^At) and perform another exponential fit 

to Cpg{At) within the plateau region in order to determine the singlet pseu- 
doscalar mass. 

We find this dual filtering approach to operate very well in the two-flavour 

case, on the SESAM QCD configurations. This is iUustratcd in Fig. (12) which 
manifests very clearly a plateau in the local mass plot of the ry'-correlator. The 
onset of plateau formation is precocious in the sense that we still maintain very 
precise signals. Therefore the statistical accuracy is sufficiently high to resolve 
the mass gap between the singlet and nonsinglet masses, whose lattice values for 
these sea quark masses are also included in the plot, being indicated with the 
symbol V. 

Let us emphasize that we consider the upshot of precocious asymptotia in 
the r/'-correlator as an essential outcome of our analysis; it corroborates that the 
nonperturbative lattice simulations take efi'cct in implicitly creating the mass 
gap between singlet and nonsinglet states, much in line with the perturbative 
reasoning from Eq. (22): indeed, the simulation data for the two-loop correlator, 
D(t), turn out to be the difference between just two exponentials rather than a 
multiexponential superposition. This will enable us to perform a rather detailed 
physics analysis in the singlet pseudoscalar channel. 

One might still worry about possible systematic errors from spectral trun- 
cation in the singlet mass plateaus contained in Fig. (12). In this context it 
is reassuring to observe nice agreement within errors with the results obtained 
from stochastic sources, as illustrated in Fig. (13). 
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Fig. 12. Mass plateau formation in the 77' (singlet) correlator after dual filtering with 
TEA, for two different values of hopping parameters (sea quark masses) versus time 
separation between source and sink, At [24]. The plot shows the respective lattice 
values for the pion (nonsinglet) mass for comparison. The effect of smearing onto the 
signal is also illustrated: SS and LL refer 'to smeared source and sink' and 'local source 
and sink', respectively. 
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Fig. 13. Comparing mass plateaus of the singlet correlator obtained with SET and 
TEA [24]. Smeared sources and sinks are used throughout. Notations as in Fig. 12. 



It is interesting to compare the plateau formation from dual filtering (Fig. (12)) 
with the state-of-the-art results from optimized smearing as obtained in the most 
recent comprehensive studies of the CP-PACS collaboration [20] (Fig. (4.3)). The 
comparison teaches us that dual filtering is indeed remarkably effective in un- 
ravelling the mass plateaus in the pseudoscalar singlet mesons. 
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Fig. 14. Statc-of-thc-art mass plateau formation from optimized smearing as obtained 
in a comprehensive study by the CP-PACS group, at their smallest lattice spacing, 
a = .1076 fm, on a 24^ x 48 lattice, for various sea quark masses. Quotation from 
Ref. [20]. 



5 Towards realistic physics results 
5.1 Unquenched two flavour world 

Given the good quality of signals within the SESAM setting, we are prepared to 
compute the entire mass trajectory of singlet pseudoscalars. In table 1 we list the 
essential run parameters of the SESAM project, which is based on vacuum field 
configurations from standard Wilson fcrmions on a 16^ x 32 lattice [8], with two 
active sea quark flavours, at one value of the coupling, /3 = 5.6. Four different 
sea quark masses have been used. The most interesting control parameter is the 
ratio between the pion and p-meson masses as determined on the lattice; it is 
quoted in the second column of table 1. 

This setting allows for a chiral mass extrapolation of the singlet pseudoscalar 
meson composed of u and d quarks, as plotted in Fig. (15). The plot contains the 



Table 1. Simulation parameters used at /3 = 5.6 and numbers of stochastic sources 
Ns (see Eq. 47 and Fig. 13). Last column: numbers of available decorrelated vacuum 
field configurations, Nconf- 





*T N, 


■^conf 


0.1560 0.834(3) 16^ 


* 32 400 


195 


0.1565 0.813(9) 16'' 


* 32 400 


195 


0.1570 0.763(6) 16^ 


* 32 400 


195 


0.1575 0.692(10) 16^ 


* 32 400 


195 



The Computation of ij and rj' in LQCD 



21 



Quadratic chiral extrapolation, 2 flavours 
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Fig. 15. Plot of versus quark mass, all given in lattice units [28]. 

mass trajectory for the singlet channel as well as the nonsinglet extrapolation, 
marked by V. We conclude that the data allow for a safe chiral extrapolation 
in both channels, though there is definite need to lower the sea quark masses 
in the simulations. The precision is sufficient to resolve the mass gap between 
singlet and nonsinglet pseudoscalars. 

As we have excluded strange quarks altogether, we can of course not expect 
to hit the experimental 77' mass: the extrapolated singlet mass of 290 MeV should 
be compared to the rj rather than to the r}' mass. It is evident, that we must 
address the issue of treating strange quarks in addtion to u and d quarks. 

Another shortcoming of SESAM is its restriction to one coupling (3 which 
does not allow for a continuum extrapolation. Very recently, the CP-PACS col- 
laboration has found in their Nf = 2 study with an improved action that this 
extrapolation tends to induce a considerable increase in the pseudoscalar singlet 
mass [20]. 

5.2 Partially quenched scenario for the three flavour world 

One of the main restrictions of todays unquenched QCD simulations is their 
limitation to two mass degenerate dynamical sea quark flavours. This presents 
a serious shortcoming when dealing with the real physics of the rj-rj' system. 
Nevertheless, the partially quenched approach [30] offers a crutch to include the 
strange quark sector into the analysis. In the partially quenched scenario, strange 
quarks might appear as valence quarks only and not contribute at all dynamically 
to the quark sea. 

For the sake of discussion, let us depart from the fully quenched situation. 
In this setting, the s-loop s-loop correlator, Dgs in Eq. (38), amounts to the 
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occurence of a double-pole: 

in momentum space. In this expression, rrig denotes the valence approximation 
mass estimate for the pscudoscalar ss meson that can be determined on the 
lattice from the connected piece of the Cgs correlator. As we have seen in section 
2.2, the double pole impedes an exponential decay in t and hence will require a 
different strategy of analysis: it can readily be rewritten as a derivative on the 
single pole expression 

2 // 2 I 2\2 ^ ™0,ss /_„^, 

mo,J{p +mj =-^^T^- (53) 

Iff tuq is a constant w.r.t. p^, Eq. (53) is readily Fourier transformed and pre- 
dicts a linear t-dependence in the ratio of disconnected to connected correlators 
to look for: 

m--=%f^=mlJ{'^rn,)xt. (54) 

Note that prior to the exploitation of this ratio relation on actual data, ground 
state filtering the Cgg-data is highly advisable, as explained in section 4.3. The 
effective strength of the double pole, mg ^g, can be extracted from Eq. (54) and 
is interpreted as the mass shift due to OZI-rule violation dynamics in channel 2, 
much in the spirit of Eq. (24): 

= + ™o,ss • (55) 

Let us turn now to the partially quenched setting. We wish to consider a situation 
with broken SU (3) flavour symmetry, where alle valence and sea quarks can differ 
in masses. 

According to Ref. [30] partial quenching of a particular quark species is 
achieved by introduction of scalar pseudoquark partners (of equal mass) that 
freeze those quark degress of freedom inside the determinant. In this context, 
one introduces a basis of states with Nq quarks and k additional pseudoquarks. 
In our situation, we have Nq = 5 (three valence and two sea quarks) and A; = 3, 
if we allow for different masses of sea and valence quarks. The correlator has to 
make reference to these degrees of freedom, In Fourier space, it is modeled 
by the momentum space ansatz [15] 



^^^•-p2+^2 (p2+^2)(p2+^2)^(p2) , i, j - ■ ■ ■ Nq + k . (56) 

The quantity Ci is equal to +1 (-1) for quark (pseudoquark = determinant eater) 
channels. The function F{p^) can be viewed as to model the implicit multiloop 
contributions from sea quarks to the disconnected pole terms: 

F{p') = l+ (57) 
, p^ + mi. 
sea quarks a 
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. 16. Our points of simulation in the hopping parameter plane of valence/sea quarks. 



rrii, rrij and rrig- denote neutral pseudoscalar nonsinglet meson masses, while mg 
refers to the effective singlet interaction that might depend on the quark species 
involved. 

F reads in the case of the SESAM simulation with degenerate dynamical u 
and d quarks 

F{p') = 1 + . (58) 

Note that the model incorporates both the fully quenched limit, where F{t) = 1, 
as well as the consistently Nj — 2 unquenched situation, where u and d valence 
and dynamical quark masses coincide, such that mi = 7712 = in this latter 
instance, the double pole term in the nonstrange sector is lifted and turned into 
a displaced single pole term, in accordance with the considerations presented in 
section 4.3. 

For practical matters, a partially quenched lattice analysis of Eq. 56 is carried 
out in two steps as follows 

1. The pseudoscalar meson masses, mi, m2, mg., are determined first by LQCD 
from connected pseudoscalar correlators with the appropriate valence quark 
settings (in the case of and mo^o-, sea and valence quark masses are 
identified) , 

2. the singlet correlator matrix contains a set of mass gap parameters, TOo,ij- 
They can be computed finally from a fit of the Fourier transform of the 
ansatz, Eq. (56), to the lattice data for C{t). 



The latter fitting can be done time-slicewise, yielding fit parameters for 'local' 
effective mass gaps which should exhibit mass plateau formation. It goes without 
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saying that throughout all the numerical analysis steps one should use the dual 
filtering method with 'synthetic' data in the sense of section 4.3. 



5.3 Partially quenched results 

In a comprehensive study [29] we analyzed some 800 independent SESAM con- 
figurations on 16^ X 32 lattices at /3 = 5.6 [31], with four different sea quark 
masses. In Fig. 16 we have plotted the set of hopping parameters of valence 
quarks in perspective with their critical values, for our different sea quark set- 
tings. The full QCD situation with 2 quark flavours in the sense of chapter 5.1 
is indicated by the dotted line, Kgea = '•^valence and the thin line connects the 
hopping parameters for strange quarks. 



m_0,nn mass plateaus, kappa_sea=0.1575, alpha=0.0 
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Fig. 17. First evidence of plateau formation in a partially quenched setting, illustrated 
on the effective mass gap, monn{t), in the disconnected two loop correlator, Dnn{t). 
The data refer to the lightest sea quark mass of SESAM, i.e. Ksea ~ .1575, for various 
isosinglet valence quark masses. 



We analyse the data according to Eq. 56. In doing so we proceed to study 
effective mass gaps as determined from the hairpin diagrams, Dij{t), by fitting, 
w.r.t. rriQijit) := /i at each value of t, the zero-momentum Fourier transforms of 
ansatz Eq. 56 to our data. In broken SU(3) with two active quark flavours, we 
introduce two isosinglet pseudoscalars, designated by the indices n (for 'light', 
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Chiral extrapolation for elements of quadratic mass matrix 
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Fig. 18. Extrapolations of the quadratic mass matrix elements of light (n) and strange 
(s) isosingiet valence quarks to QCD with chiral sea quarks (Eq. 59), at a = 0. 



nonstrange) and s (for strange). In this manner, we can study three types of 
hairpin diagrams, -D„„ , Dns , D^s, with the strange quark mass determined 
on the lattice from the kaon mass (for k locations, see Fig. 16). As a result, 
we find satisfying plateau formations in mQij(t) as illustrated in Fig. 17. In 
fact, the numerical quality of our signals comfortably allows for the consecutive 
extrapolations (i) to light valence and then (ii) to light sea quark masses; the 
chiral sea quark extrapolation is exhibited in Fig. 18. 

From the consecutive chiral extrapolations, we obtain for the quadratic mass 
matrix in the quark-flavour basis 



M 



2 



(59) 



where m^„ :— A/^„ + mQ„„. With the p- meson mass scale [8] we find 

/ (306 ±91)2 72(201 ±34)2. 
x/2(201 ± 34)2 (ggo ± 15)2 I MeV 



Diagonalization of A4 renders 



Mr, = 292 ± 31 MeV, Mr,> = 686 ± 31 MeV. 



(60) 
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Fig. 19. Minimizing the valence quark dependence of the mass gap by variation of a. 
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These numbers from Nf = 2 look promising: the individual mass values - while 
being low w.r.t. the real three flavour world - yield a mass splitting that compares 
very well to phenomenology. 



5.4 Data in accord to xPT? 

So far we have neglected that the partially quenched scenario has been devised 
in the framework of chiral perturbation theory (xPT) where we expect to deal 
with a limited number of effective couplings. In fact, Bernard and Golterman [30] 
showed that the effective chiral symmetry violation (due to the chiral anomaly) 
is most simply described by the following contribution to the chiral Lagrangian: 



Co - 



2 



l\2 



(61) 



i.e. in terms of just two constants, /i (mass gap in the chiral limit) and a (in- 
teraction parameter) only! If the tree approximation of x PT to their partially 
quenching scenario holds, the simulation data should be fitted in terms of these 
two constants, irrespective of the (supposedly light!!) valence quark masses cho- 
sen. We will now address the question whether this is indeed the case. To this 
end we have to first rewrite the correlators, Eq. 56, by including the parame- 
ter a. Actually, the Lagrangian translates into two-point correlation functions 
(propagators) for the neutral pseudoscalar quark bilinears 



'^ri{x) = gi(a;)75(ji(x) , 



(62) 
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Fig. 20. The analogue of Fig. 17: Universal plateau formation, at the lightest sea quark 
mass, when adjusting a to the value 0.0328. 



In momentum space the disconnected part of these ps-correlators again has 
a compact form [32] with the infamous double pole structure 

The prefactor, TV = 1/(1 + A^/Q!) and the singlet mass, M^, = {Mj + Nfn^)/{1 + 
Nfa) are seen to carry explicit a-dependencies. Rememeber that the octet 
masses {Mii,Mjj, and Mj) from valence and sea quarks, respectively, are readily 
computed from connected diagrams on the lattice. Our present concern is the 
determination of the mass gap, and a from our lattice data. 

We remark that the outcome of the analysis at this stage seems not to be 
in accord with the tree approximation to the Lagrangian, Eq. 62, as the latter 
predicts the mass plateaus to be independent of the valence quark mass, my, 
contrary to the apparant findings in Fig. 17. Therefore, we repeated the above 
mass plateau study on a set of nonvanishing a-parameters, in an attempt to 
verify the compliance of our data with such independence: inspecting the plot 
of rrionn vs. niy (Fig. 19) it becomes obvious (i) that m^nni^v) shows a simple. 
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linear behaviour and (ii) that a can indeed be adjusted to produce a zero slope 

of monnimy) 

In this way it is straightforward to find an optimal value, aopt, that elimi- 
nates (for a given sea quark mass) the m^-dependence from monn, resulting in a 
universal plateau level. This collapse of data into universal plateaus (of heights 
jl) is exemplified in Fig. 20. The comparision with the situation encountered 
with a = (Fig. 17) provides clear evidence for our sensitivity in determining 

O^opt • 

After chiral sea quark extrapolation wc arrive at a first estimate of a: 

a = 0.028 ± 0.013 ^ = 203 ± 34MeV . (64) 

Note that the numerical value for the mass gap at the chiral point is robust w.r.t. 
the two different approaches presented here. One might phrase it like this: the 
allowance for a nonvanishing a parameter is tantamount to the anticipation of 
the respective chiral valence quark limit in the mass gap! 

6 Conclusion 

The application of spectral techniques provides unprcccdcntcdly high accuracy in 
the study of hairpin diagrams. Spectral methods thus provide access to detailed 
studies of the flavour singlet mesons. In this work, this has been demonstrated 
in two flavour QCD as well as in the A^^ = 2 partially quenched scenario, with 
standard Wilson fermions and in the regime of medium mass sea quarks (see 
Fig. 16. Since spectral methods are at their best in the truly chiral regime, this 
looks very promising in regard to future simulations of flavour singlet mesons 
with realistically light sea quarks, that we expect to come up within the overlap 
fermion formulation on the lattice. 
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